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We present a theoretical method to study driven-dissipative correlated systems on lattices with 
two spatial dimensions (2D). The steady-state density-matrix of the lattice is obtained by solving the 
master equation in a corner of the Hilbert space. The states spanning the corner space are determined 
through an iterative procedure, using eigenvectors of the density-matrix of smaller lattice systems, 
merging in real space two lattices at each iteration and selecting M pairs of states by maximizing 
their joint probability. Accuracy of the results is then improved by increasing M, the number of 
states of the corner space, until convergence is reached. We demonstrate the efficiency of such an 
approach by applying it to the driven-dissipative 2D Bose-Hubbard model, describing, e.g., lattices 
of coupled cavities with quantum optical nonlinearities. 

PACS numbers: 


Simulating large quantum systems is a challenging task 
because their complexity grows exponentially with their 
size. Indeed, the dimension of the Hilbert space for a 
multipartite system consisting of m subsystems, each of 
them described by a space of dimension N, is N^. Fur¬ 
thermore, for open systems the physics can no longer be 
described only in terms of the eigenstates of an Hamil¬ 
tonian, requiring instead the knowledge of the density- 
matrix. In this case, the number of variables to be de¬ 
termined scales as namely the square of the size of 

the Hilbert space. 

In the last decades, several methods have been pro¬ 
posed to reduce the complexity of this problem. The first 
attempt in this direction is the renormalization group 
technique, proposed by Wilson [1] and successfully ap¬ 
plied to the Kondo problem. Numerical implementations 
of this approach are based on the solution of a system 
with a smaller Hilbert space, where only the relevant 
physical states with the lowest energies are retained. Ide¬ 
ally, this procedure can be iterated by arbitrarily growing 
the size of a block system step by step, for instance by 
doubling the size of the block at each iteration. How¬ 
ever, such a numerical implementation of the real-space 
renormalization group can yield inaccurate results for 
the system ground state, because the boundary condi¬ 
tions imposed while solving the smaller system might be 
inappropriate to describe the doubled one [2]. In the 
case of one-dimensional systems, a powerful method is 
represented by the density-matrix renormalization group 
(DMRG) [3], which is based on the selection of the most 
probable states of the reduced density-matrix of a block, 
obtained by determining the ground state of the Hamilto¬ 
nian of a larger section of the lattice. The generalization 
to two spatial dimensions is challenging and currently un¬ 
der intense study [4, 5]: one approach exploits the artifi¬ 
cial description in terms of one-dimensional systems with 
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long-range interactions [6] , while another is based on the 
generalization of matrix product states [7] to projected 
entangled-pair states [8] . These theoretical methods have 
been extended to ID dissipative lattice systems by intro¬ 
ducing matrix product density operator algorithms [9], 
time-dependent DMRG [10] and a superoperator renor¬ 
malization technique [11]. These approaches are aimed 
at solving the master equation governing the dynamics 
of the density-matrix of the lattice. 

Among driven-dissipative systems, lattices of cou¬ 
pled cavity resonators with quantum optical nonlinear¬ 
ities [12-14] are attracting a considerable interest, e.g. 
for the realization of non-equilibrum strongly correlated 
photonic phases [15]. In particular, the spectacular rise 
of circuit QED resonators with superconducting Joseph- 
son quantum circuits is very promising in this respect 
both for the realization of strong correlations and for 
their control [16, 17]. So far, several studies have been 
devoted to non-equilibrium mean-field-like theories [18- 
23], based on a Gutzwiller factorization of the density- 
matrix. Numerical methods beyond mean-field for such 
systems so far rely on a direct integration of the density- 
matrix for small size systems [24, 25] or applications of the 
matrix product operator techniques mentioned above to 
one-dimensional cavity arrays [26, 27]. 

In this letter, we present a theoretical method to ex¬ 
plore the physics of driven-dissipative correlated quan¬ 
tum systems with two spatial dimensions. A corner of the 
Hilbert space for a lattice system is selected using eigen¬ 
vectors of the density-matrix solving the master equa¬ 
tion for smaller clusters. At each step, two sublattices 
are merged and M pairs of states are selected to con¬ 
struct a corner basis by maximizing their joint proba¬ 
bility. The degree of accuracy can be controlled by en¬ 
larging the number of states of the corner space, until 
convergence is obtained. The method is applied to the 
driven-dissipative 2D Bose-Hubbard model, which de¬ 
scribes, e.g., two-dimensional arrays of coupled cavities 
with quantum nonlinearities. 

The general problem we aim to solve is the Lindblad 
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master equation [28] for the density matrix p of a driven- 
dissipative manybody quantum system, 


dp i 


^\c,0,-l{c'Ap + pCjC, 


where H is the Hamiltonian of the lattice system and Cj 
are operators describing the relaxation of the system due 
to the interaction with an external bath. We will consider 
a zero-temperature reservoir for simplicity, although the 
case of finite temperature can be treated without major 
complications. To give an example, in the case of a lattice 
of optical cavities the jump operators are Cj = 

where bj is the photon annihilation operator in the j-th 
cavity and 7 j is the corresponding dissipation rate. 

The corner space renormalization algorithm we intro¬ 
duce here is based on the following steps: i) determine the 
steady-state density-matrix for small lattices, for which 
a direct, brute-force integration of the master equation is 
possible; ii) merge spatially two pre-determined lattices 
and select the M most probable product states spanning 
the so-called corner space; hi) determine the steady-state 
solution of the density-matrix in the corner space; iv) in¬ 
crease the dimension M of the corner until convergence of 
the observables is achieved; v) in order to create a larger 
lattice, go back to step ii). 


Here, we describe in detail the crucial steps hi) and 
iv), i.e., the selection of the corner space. As sketched 
in Fig. 1, let us suppose that we know the solution for 
the steady-state density matrices and for two 
lattices A and B. If we want to consider a lattice ob¬ 
tained by merging spatially the two lattices A and B, 
the corresponding Hilbert space is H(aub) = ^ 

where Ha and Hb are the Hilbert spaces of A and 
B. Each density-matrix operator can be diagonalized 
as p^-^^ = where the states 

form an orthonormal basis for Ha and are the cor¬ 
responding probabilities. Analogous notations hold for 
the system B. Each ket \(t)i^^) represents a manybody 
state which can have strong correlations within the sys¬ 
tem A. To select a small ‘corner’ C{M) of the larger 
space Haub^ we will consider the subspace spanned by 
the M most probable states of the form 
ranked according to the joint probability Let 

us call the most probable product state, i.e. 

such that ^ cvcry value of r and 

r'. We will call the second most probable 

product state and so on so forth. In other words, we 
will consider the subspace generated by the orthonormal 
basis )>•••> where 
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FIG. 1: Sketch of the corner space renormalization method. 


M 

n 

mb)) 

9<i,l> 

20 

0.09443 

0.2772 

1.029 

50 

0.09469 

0.2770 

0.9693 

100 

0.09513 

0.2768 

0.9652 

200 

0.09541 

0.2767 

1.061 

400 

0.09544 

0.2767 

1.058 

800 

0.09549(3) 0.27671(5) 

1.0644(1) 

1600 

0.09547(3) 0.27672(6) 

1.0643(1) 

65536 

0.0954(1) 

0.2764(2) 

1.0643(3) 


TABLE I: Corner method results for the driven-dissipative 
Bose-Hubbard model with periodic boundary conditions and 
the following parameters: 4x4 square lattice {z = 4), 
U = +00 {Nmax = 1, hard-core bosons), J /7 = 1, F/j = 2, 
Au/j = 5. The numbers in parenthesis indicate the statis¬ 
tical errors on the last significative digit due to finite Monte 
Carlo sampling when applied. In this example, the dimen¬ 
sion of the full Hilbert space is 2^® = 65536. The case of 
65536 states has been solved by an independent Monte Carlo 
wavefunction code using a Fock basis for the entire space and 
sparse matrix calculations. 


M most probable pairs of states^. Note that a generic 
state belonging to the corner space, namely of the form 
1^) = Xlfli describe strong correla¬ 

tions and quantum entanglement between systems A and 
B while keeping correlations within A and B. We empha¬ 
size that by increasing arbitrarily the number M of states 
in the corner space, the method becomes exact, because 
the considered basis spans the entire Hilbert space. Of 
course, the method is useful only when the number of 


^ Note that can be smaller than or can be smaller 

than ^ , but always pif^p^J^ > ■ 
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M 

n 

mb}) 

92 

s<j,i> 

20 

0.0902 

0.1967 

1.646 

1.28 

50 

0.1006 

0.1907 

1.513 

1.34 

100 

0.1044 

0.1886 

1.454 

1.26 

200 

0.0968 

0.1922 

1.324 

1.51 

400 

0.1006 

0.1905 

1.291 

1.51 

800 

0.1009(2) 

0.1903(3) 

1.242(3) 

1.57(2) 

1600 

0.1014(2) 

0.1896(2) 

1.226(3) 

1.58(2) 

3200 

0 . 1002 ( 2 ) 

0.1897(2) 

1.185(2) 

1.63(2) 

6400 

0.0994(2) 

0.1899(2) 

1.179(3) 

1.63(1) 


TABLE II: Parameters: 4x4 square lattice with periodic 
boundary conditions, V— 20, Jj'y — 3, Fj^ — 2, = 

5. A number Nmax = 3 of bosons per site has been considered. 
In this case, the dimension of the full Hilbert space is 4^® ^ 
4.3 • 10^ 


states M required to reach convergence is small enough 
to be treated numerically. This ultimately depends on 
the degree of correlation of the considered system. 

Concerning step iii), namely the determination of the 
steady-state density-matrix, it is worth pointing out that 
when the number M of states in the corner space is small 
enough, the master equation can be solved by direct nu¬ 
merical integration in time (M typically up to a few hun¬ 
dreds). For larger values of M (up to a number of the 
order of 10^ depending on the sparsity of the Hamilto¬ 
nian matrices), a more efficient method is based on a 
stochastic technique [29], the so-called Monte Carlo wave- 
function algorithm [30-32] . Such algorithm computes the 
density-matrix of the system by averaging over quantum 
trajectories of the wavefunction in presence of random 
quantum jumps. 

As a first illustration of the corner space renormaliza¬ 
tion method, we show results for the driven-dissipative 
Bose-Hubbard model in 2 D square lattices. The cor¬ 
responding Hamiltonian {h = 1) in the frame rotating 
at the pump frequency and in the case of homogeneous 
pumping reads: 



t(1/y) 


FIG. 2 : Evolution of n and g 2 versus time t (units of I/ 7 ) 
for the driven-dissipative Bose-Hubbard model with periodic 
boundary conditions on lattices of various size for the follow¬ 
ing parameters: Uj^ — 20, J /7 = 3, F /7 = 2 , Alu/j = 5 . 
Solid lines represents evolutions performed by direct integra¬ 
tion of the master equation, while points depict Monte Carlo 
wavefunction calculations. When error bars are not shown, 
the statistical error is smaller than the point size. The black- 
dotted lines represent the mean-field values. The initial con¬ 
ditions are explained in the text. 



State Rank r State Rank r 


H = Y,{-Au:^bJ + ^^^bJbJ+Fi^+bJ))F 

j <j,l> 

where Auj = uOp — uOc is the detuning between the pump 
and the bare boson frequency, U is the on-site boson- 
boson interaction and F is the pump field. J is the hop¬ 
ping coupling, z is the coordination number and 
denotes the sum over all the couples of nearest neighbors. 
For simplicity, we have fixed the phase of the pump in 
such a way that F is real. Finally, each site is subject to 
losses with a dissipation rate 7 . 

In the following, we will consider the case of periodic 
boundary conditions. In Table I we show results for a 
4x4 square lattice of hard-core bosons (U = +oc), i.e. 
for which the maximum number of photons per site to 


FIG. 3: Probabilities pr (top panels, logarithmic scale) and 
expectation value of the total boson population (ritot) = 
^j{nj) for the orthonormal eigenvectors |Tr) of the steady- 
state density-matrix (p = pr|4^r)(4^r| and pr > Pr+i). 

The state rank r is in logarithmic scale. Lattice size: 6x3. 
Driving parameters: F /7 = 2, Alu/j = 5. Left: Uj^ — 2^ 
and J /7 = 3. Right: hardcore bosons with J /7 = 1. 


be considered is Nmax — 1- These results have been ob¬ 
tained starting from a 2 x 2 lattice for which a brute-force 
determination of the steady-state solution of the master 
equation is possible. Merging two 2x2 lattices, we get 
results for a 4 X 2 lattice and, repeating the doubling 
procedure, for the 4x4 case. The dimension of the full 
Hilbert space for hard-core bosons on a 4 x 4 lattice is 
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2 ^^ = 65536. Although very heavy, the master equation 
resolution in the full Hilbert space has been performed 
by an independent Monte Carlo wavefunction code using 
a Fock basis for the entire space and sparse matrix cal¬ 
culations. This way, we have been able to benchmark 
the results obtained with a small corner to the exact 
results. The table reports results for the boson pop¬ 
ulation per site n = (5^6^), the real part {^{{bj))} of 
the bosonic coherence and the nearest-neighbor correla¬ 
tion por hard-core bosons, the on-site 

rijUi ’ 


second-order correlation tunction g 2 = ^ ^ 2 —- is triv¬ 

ially equal to 0 since two bosons are not allowed to be 


on the same site. Note that = 1 for a factorized 

Gutzwiller-like density-matrix pQ = pj where pj is 
the reduced density-matrix of the j-th site. The mean- 
field approach is equivalent to taking a self-consistent 
Gutzwiller density-matrix, assuming all the sites identi- 
cal. Hence the difference {g)^j — 1) quantifies the de¬ 

gree of correlations beyond mean-field between nearest 
neighbors. Remarkably, for the parameters given in the 
caption of Table I, we get a very accurate result for a 
small number M = 200 (negligibible error for n, 0.1% 


for the bosonic coherence and 0.3% for In Ta¬ 


ble H, we show results for soft-core bosons with a larger 
hopping coupling {U = 20, J /7 = 3) and a cut-off 
number N^ax = 3 of bosons per site (we have verified 
that this is the cut-off number per site required to get 
convergence). A cut-off Nmax = 3 for a 4 x 4 lattice im¬ 
plies a Hilbert space dimension equal to 4}^ 4.3 • 10^. 

As shown by the convergence progression in Table H, re¬ 
sults with deviations below 1 % are reached already for 
a corner space dimension M = 3200, hence six orders 
of magnitude smaller than the full Hilbert space for a 
system exhibiting large correlations — 1 = 0.63). 


An example of the temporal dynamics leading to 
steady-state solutions is reported in Fig. 2, plotting n 
and g 2 for different lattice sizes. The corner method re¬ 
sults are compared with the non-equilibrium mean-field 
approach used in Refs. [ 21 , 23], based on the exact an¬ 
alytical solution of the master equation for the one-site 
problem [33]. The initial condition for the density-matrix 
dynamics for the 2 x 2 lattice is the mean-field solution. 
After a transient, a steady-state solution is obtained. The 
initial condition for the 4x2 lattice is constructed from 
the steady-state solution of the 2 x 2 lattice and so-on 
so forth. We have also merged 3 x 1 clusters to get the 
3x3 lattice and then the 6x3 case by doubling. We 
see that the steady-state observables for the 3 x 3, 4 x 4 
and 6x3 lattices with periodic boundary conditions tend 
to converge to the same value, so the results are already 
approaching those for a lattice with an infinite number 
of sites. The finite spatial range of the correlations of 
the driven-dissipative system is responsible for such rel¬ 
atively quick convergence. For the parameters in Fig. 
2 , the deviations from the mean-field theory are around 



Mean-field 

Corner method 

Uh 

n 



n 


9(3,1) 

00 

0.0953 

0 

g X 4 ( 1600 ) 

8 X 

0.09527(2) 

0.0948(2) 

0 

0 

1.0436(3) 

1.0237(6) 

20 

0.125 

0.836 

4 X 

6 X 

0.1281(4) 

0.1282(9) 

0.859(4) 

0.858(9) 

1.172(5) 

1.173(4) 

20 * 

0.0768 0.8879 

4 X 4 ^®^°°) 

6 X 

0.0994(2) 

0.0992(1) 

1.179(3) 

1.202(4) 

1.63(1) 

1.65(1) 

10 

0.9587 0.6088 

4 X 

3 X 3(«“°) 

0.9275(8) 

0.9281(9) 

0.631(1) 

0.617(1) 

1.0127(8) 
1.0069 ( 6 ) 

1 

0.1156 

1.265 

16 X 8^®°°) 

0.1156 

1.259 

0.9897 

0.5 

0.1126 

1.112 

16 X 16^'‘°°) 

0.1126 

1.1105 

0.9941 


TABLE III: Steady-state expectation values for lattices (peri¬ 
odic boundary conditions) with different sizes, calculated via 
the Gutzwiller mean-field theory and the corner space renor¬ 
malization method. M is the dimension of the corner space. 
Parameters: J\ (except the third line with the * sign, 
obtained with J/7 = 3), F/7 = 2 and Acj/7 = 5. The max¬ 
imum number of bosons per site is Nmax = 1 for hardcore 
bosons, Nmax = 3 for F/7 = 20 , Nmax = 5 for F/7 = 10 , 
Nmax = 4 for F/7 = 1 and 0.5. 


20 % for n and ^ 2 - Since the driving is homogeneous and 
the considered boundary conditions are periodic, short¬ 
comings due to conflicting boundary conditions in the 
doubling procedure do not apply here. 

It is insightful to look at the diagonal decompo¬ 
sition of the calculated density-matrix, namely p = 
Pr|Tr)(T^| where Pr > Pr+i- In Fig. 3, we show 
an example of the probability distribution pr (top panels, 
logarithmic scale) and the expectation value of the total 
number of bosons in the lattice (bottom panels) versus 
the state rank r for a 6 x 3 lattice of soft-core bosons 
with U = 2 O 7 (left panels) and hard-core bosons (right 
panels). In both cases, the probability drops sharply by 
several orders of magnitudes when the rank r is large 
enough, confirming the achieved convergence of the cor¬ 
ner space dimension. In the hard-core boson case, a 
rather well definite shell structure is apparent. The first 
state (r = 1 ), which captures a large part of probability, 
is followed by shells of states having close probabilities 
and densities. In the case of a homogeneous system, a fac¬ 
torized Gutzwiller density-matrix with each site having 
the same reduced-density matrix leads to a shell struc¬ 
ture with exactly flat plateaux structures due to sym¬ 
metry reasons. In fact, a permutation of the role of the 
different sites does not change the probability of a state 
and observables like which is a sum over all the 

sites. In the right panel of Fig. 3 (hard-core boson case), 
the situation is qualitatively close to the Gutzwiller case, 
even though the plateaux are not exactly flat. In the 
case of soft-core bosons in the left panel of Fig. 3, a first 
shell is clearly visible, while higher shells merge into a 
continuous curve where the different quantities increase 
gradually, denoting a large degree of correlations (indeed 
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— 1 ^ 0.6 in the case considered). 

In Table III, we summarize illustrative results for dif¬ 
ferent lattice sizes with periodic boundary conditions and 
compare them to the Gutzwiller mean-field solutions [21, 
23], using the same excitation parameters (T/y = 2 and 
Acj/y = 5). The convergence of the results with increas¬ 
ing corner dimension M has been checked as well as the 
required maximum number N^ax of bosons per site. It is 
apparent that in the considered case the mean-field the¬ 
ory gives rather good results for hard-core bosons and a 
large 8x8 lattice, as quantified by a — 1 ^ 0.02. Sig¬ 
nificant deviations are instead present when the on-site 
interaction U is competing with the hopping coupling J 
(the cases with I//7 = 20 and J/7 = 1 and 3 in Table 
III). The value for I//7 = 10 and J/7 = 1 is close to a 
two-photon resonance [23] and indeed the the population 
of bosons per site is much higher (close to one boson per 
site) with the on-site g 2 correlation function quite close 
to 0.5. For I//7 = 0.5, it is possible to simulate very 
large lattices (a 16 x 16 lattice is reported) with a very 
small number of states (M = 400). 

In conclusion, we have presented a theoretical method 


for driven-dissipative 2D correlated lattice systems. The 
proposed numerical algorithm follows a hybrid real-space 
renormalization group approach where the states are se¬ 
lected on the basis of joint probabilities. We have suc¬ 
cessfully demonstrated the efficiency of such a method 
by applying it to the driven-dissipative Bose-Hubbard 
model on 2D square lattices. Unlike mean-field theories, 
where the decoupling approximation is not controlled, 
the present numerical method allows us to get results 
with controllable accuracy. The method has therefore 
the potential to become a precious tool to benchmark 
analytical theories and study strongly correlated open 
systems with more than one spatial dimension. Future 
studies will explore the physics of 2D arrays of nonlinear 
cavities with complex elementary cells (including disor¬ 
der), geometric and spin frustration as well as the role of 
artificial gauge fields in extended lattices. 

C. C. acknowledges support from ERG (via the Gon- 
solidator Grant ’GORPHO’ No. 616233), from ANR 
(projects QPOL and QUANDYDE) and from Institut 
Universitaire de Erance. 
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